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Abstract 

In the study of open quantum systems, the polaron transformation has recently attracted a 
renewed interest as it offers the possibility to explore the strong system-bath coupling regime. 
Despite this interest, a clear and unambiguous analysis of the regimes of validity of the polaron 
transformation is still lacking. Here we provide such a benchmark, comparing second order pertur- 
bation theory results in the original untransformed frame, the polaron frame and the variational 
extension with numerically exact path integral calculations of the equilibrium reduced density ma- 
trix. Equilibrium quantities allow a direct comparison of the three methods without invoking any 
further approximations as is usually required in deriving master equations. It is found that the 
second order results in the original frame are accurate for weak system-bath coupling, the full 
polaron results are accurate in the opposite regime of strong coupling, and the variational method 
is capable of interpolating between these two extremes. As the bath becomes more non-Markovian 
(slow bath), all three approaches become less accurate. 

PACS numbers: 03.65.Yz, 71.38.-k 
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I. INTRODUCTION 



In many open quantum systems, the coupling between the system and the bath can be 
considered as a small parameter. In this case the application of second order perturbation 
leads to a master equation of Redfield or Lindblad type^. Their numerical implementation is 
straightforward and not computationally expensive. However, for many physical systems of 
current interest it has been shown that the weak coupling approximation is not justified. One 
example is the energy transfer process in photosynthetic complexes where the magnitude of 
the system-bath coupling is comparable to the electronic couplings^i"-. There are only a few 
non-pert urbative techniques to obtain the numerically exact dynamics; examples include the 
hierarchy master equation^''', the quasi- adiabatic propagator path integral (QUAPI)^'^, and 
the multiconfiguration time-dependent Hartree (MCTDH)^"— approach. However, these 
methods are computationally demanding and also not trivial to implement. 

Recently, a polaron transformed second order master equation^^"^'' and its variational 
ioim^^ have been derived to study the dynamics of open quantum systems at strong 
coupling. This approach transforms the total Hamiltonian into the polaron frame such that 
the system Hamiltonian is dressed by a polaron. The master equation is then obtained 
by applying perturbation theory to the transformed system-bath interaction term. This 
approach extends the regime of validity of the master equation to stronger system-bath 
coupling, provided that the electronic couplings (or tunneling matrix elements) are small 
compared to the typical bath frequency. When this condition is not fulfilled the polaron is 
too sluggish to accurately follow the system motion and the polaron transformation may 
perform worse than the standard master equation approach. 

In order to partially overcome this difficulty, the variational method has been developed 
as a generalization of the polaron transformation^^'^^. Instead of performing the full trans- 
formation, the variational polaron approach seeks for an optimal amount of transformation, 
depending on the properties of the bath. Thus it is able to interpolate between the strong 
and weak coupling regimes and to capture the correct behavior over a much broader range of 
parameters. Both the polaron and variational master equations have the attractive feature 
of being computationally economic (they have the same computational complexity as the 
Redfield equation) and are therefore suitable for studying large systems. 

However, a thorough assessment of the accuracy of second order perturbation theory in 
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the polaron and variational polaron frames is still lacking. It is not exactly clear how the 
accuracy depends on the properties of the bath, namely the bath relaxation time and the 
coupling strength. One of the main goals of this work is to provide such a benchmark. In- 
stead of studying the dynamics, here we focus on the equilibrium density matrix. Focusing 
on this quantity offers two key advantages. Firstly, in the equilibrium case the second order 
perturbation is the only approximation involved. In the derivation of second order master 
equations, additional approximations generally must be invoked, such as factorized initial 
conditions, the Born-Markov approximation, the rotating wave approximation, etc. These 
additional restrictions prevent a clear assessment of the isolated role of second order pertur- 
bation theory and the merits of the polaron transformation. Thus studying the equilibrium 
density matrix offers a direct comparison of the various perturbation methods. A second 
advantage of studying the equilibrium state is that it is much easier to obtain numerically 
exact results. Therefore we are able to systematically explore a large range of the parameter 
space that is often not possible with other exact treatments of the dynamics. 

In the next section, the details of the spin-boson model used in the remainder of the 
text are outlined. Following this, the polaron transformation and its variational extension 
are applied to the Hamiltonian in Sec. IIIBi In Sec. IlICi the second-order corrections to 
the equilibrium reduced density matrix are derived in the original, polaron and variational 
polaron frames. In the ensuing section, results for the various perturbation theories are 
compared with exact numerical results from path integral calculations over a broad range 
of the parameter space. It is found that the second order results in the original frame are 
accurate for small system-bath coupling, the full polaron results are accurate in the opposite 
regime of strong coupling, and the variational method is capable of interpolating between 
these two extremes. All three approaches become less accurate for slow baths. 

II. THEORY 

A. Spin-Boson Model 

The spin-boson model is a paradigm for the study of quantum dissipative systems. It 
has been used to investigate the energy transfer in light harvesting systems^^*^, the problem 
of decoherence in quantum optics^^, tunneling phenomena in condensed medis^^^ and 
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quantum phase transitions^IiSS. The spin-boson model consists of a two-level system coupled 
to a bath of harmonic oscillators. Its Hamiltonian can be written as (we set h = 1) 

k k 

where gi {i = x^y^ z) are the usual Pauli matrices, e is the energy splitting between the two 
levels, and A is the tunneling matrix element. The bath is modeled as a set of harmonic 
oscillators labeled by their frequencies, ujk^ and couplings to the two-level system denoted 
by gk- 

The properties of the harmonic bath are completely determined by the spectral density, 
J{(jj) = TT^^^k 9k^i^ ~ ^k)' Throughout the paper, we use a super-ohmic spectral density 
with an exponential cut-off, 

^H = ^7^e--/'*^=, (2) 

where 7 is the system-bath coupling strength and has the dimension of frequency. The cut- 
off frequency is denoted by cJc, and its reciprocal governs the relaxation time of the bath, 
r oc — . 



B. Polaron and Variational Polaron Transformation 

The polaron transformation is generated by the unitary operator 

{7 = expk5]^(6l-y, (3) 

which displaces the bath oscillators in the positive or negative direction depending on the 
state of the two-level system. The parameter fk determines the magnitude of the displace- 
ment for each mode. Setting = gk corresponds to the full polaron transformation whereas 
fk = corresponds to no transformation. The variational method allows us to determine 
an optimal value of fk that lies in between these two limits, < fk < gk^ making the 
transformation valid over a wider range of parameters. 

Applying the transformation to the total Hamiltonian, we have 

Htot = UH'totU\ (4) 
= i^o + Hi, 



where the total free Hamiltonian is Hq = Hs + Hb- The transformed system Hamiltonian 
is given by 

and the bath Hamiltonian remains unaffected, Hb = X^/e ^/c^I^/c- The transformed interac- 
tion Hamiltonian becomes 

Hi = a^Vx + ^yVy + a^V;, (5) 

where 

K = ^(^2 _25)^ (g) 

Vy = J^iD'--Dl), (7) 

v; = + (8) 

k 

and D± is the product of displacement operators, D± = H/c^^"''*^^* The tunnehng rate 
is renormahzed by the expectation value of the bath displacement operators, = BA, 
where 

trB[I)^e-^^«] 
trB[e-/3^B] 

/I 



^ = 01),.= '-^^^^^^, (9) 



exp 



2j2^coth{/3uk/2) . (10) 



Note that the interaction term is constructed such that {Hi)ho = 0. 

Following Silbey and Harris^SiSi^ we use the Bogoliubov variational theorem to determine 
the optimal values for the set {fk}- We first compute the Bogoliubov- Feynman upper bound 
on the free energy, Ab 

Ab = -^Intvs+B [e-^^°] + {Hi)h,- (H) 

Since (Hj) = by construction, the upper bound is solely determined by the free Hamilto- 
nian. The variational theorem states that Ab > A where A is the true free energy. Therefore, 
we want to make this bound as small as possible by minimizing Ab with respect to {//e}, 
i.e. = 0. The minimization condition leads to 

fk = gkFM, (12) 

r 1-1 
Fiuk) = 1 + — coth(/3a;fc/2)tanh(/3?7/2) , (13) 



where r] = + A^. 

In the continuum hmit, the renormahzation constant can be written as 

du) J{uj) 



B 



exp 



-F(cj)2coth(/3cj/2) 



(14) 



Since F{uj) is also a function of 5, the above equation must be solved self-consistently. 



C. Second Order Perturbation Theory 

This section is dedicated to finding the second order correction to the equilibrium state 
of the system. The exact equilibrium reduced density matrix can be formally written as 

Expanding the operator e~^^*°* up to second order in if/, we have 



—PHtot -/SHo 
e e 



1 _ /" dP'e^'^'^Hie-^'^" + [ d/3' [ (^/3"e^'^°if,e-(^'-'^")^°i^,e-^"^° 
Jo Jo Jo 



(16) 

The above expansion is similar to the Dyson expansion, with /3 treated as imaginary time. 

Since {Hj)^^ = 0, the leading order correction to ps is of second order in Hj. Inserting 
the above expression into Eq. ( ITSl) and keeping terms up to the second order in Hj^ the 
system equilibrium state can be approximated 

PS ^ P?^+/3?) + ...; (17) 
(2) _ _ ^-ms 



where 



A = [ dp' [ dp''Y,CU^^'-n><e^''-'^'''^ne^'''-''^^^ 
■^0 Jo 

Zf^=tvs[e-^H zP=tTs[A], (18) 



and Cnm{^) is the bath correlation function in imaginary time, 



^ , . trsfe-^^Ke-^^^y^e-^^^] 

Cnm(r) = - . (19) 

trsfe-^^s] 



The non- vanishing bath correlation functions are 



C.,{t) = -C,.(t) (20) 

= r - (21) 

^Jo TV u ^ ^ sinh(/3w/2) ' ^ ^ 

C..{r) = ^(e*(-)+e-^(^)-2), (22) 



C'^.l^) = ^(e^(^)-e-^(^)), (23) 



OO 



dio ^, cosh (i(/3 - 2t)w) 



C„(r) =/-/(.)[! -F(.)P—^;5;^^^, (24) 



where 



,M^4r^£<^nc.)-'°'\^!'f (25) 

^ Jo ^ ^ sinh(/3cj/2) ^ ^ 

It is useful at this point to analyze the behavior of the perturbation theory at strong 
coupling in the polaron frame. As seen from Eq. ( IT41) . when 7 ^ cx) then B ^ and the 
system becomes incoherent since the coherent tunneling element vanishes. At the same time, 
F{(jj) ^ 1 as S ^ so that all of the above correlation functions vanish, and hence also 
the second order correction to ps- Therefore in this limit, the equilibrium density matrix is 
only determined by the energy splitting of the two levels, ps oc exp(— 

The full polaron result can be conveniently obtained by setting F{(jj) = 1; the only non- 
vanishing correlation functions in this case are Cxx and Cyy. The opposite limit of F{(jj) = 
corresponds to performing no transformation and Czz is the only non-zero correlation 
function. For comparison, below we will also include the results from these two limiting 
cases. 



III. RESULTS AND DISCUSSIONS 



In this section, we compare the results from second order perturbation theory (2nd-PT) 
in the original [F{(jj) = 0], the polaron [F{(jj) = 1], and the variational polaron [F{(jj) as in 
Eq. f fT3l )] frames with those from numerically exact imaginary time path integral calculations. 
We compute the expectation value, (a^), since it is not affected by the transformation, 
{WazU) = {cFz)' Therefore this quantity allows us to make a direct comparison between 
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the path integral results, which provide the density matrix in the original frame, and the 
(variational) polaron results. Results from the transformed zeroth order density matrix, 
which depends only on the renormalized system Hamiltonian Hs^ are also included. 

We first calculate (a^) as a function of the dissipation strength for fast, slow and 
adiabatic baths, assessing the accuracy of 2nd-PT for different bath cut-off frequencies. 
We then conclude this section by presenting phase diagrams of the relative errors of the 
various methods as functions of the dissipation strength and the bath cut-off frequency. 
This allows us to establish the regimes of validity of each approach across the entire range 
of bath parameters. Throughout the paper, we set e = 1 and /3 = 1. 

A. Fast Bath, (jJc> A 

The value of (a^) is plotted as a function of the dissipation strength, 7, in Fig. [T]for a 
fast bath, ujc > A. Firstly, it can be seen that the result from the usual 2nd-PT in the 
original frame (dashed line) is linearly dependent on 7. While this approach is accurate at 
small 7, it quickly degrades as the coupling increases. On the other hand, the results from 
2nd-PT in the polaron (empty circles) and the variational polaron (solid line) frames are 
in excellent agreement with the exact path integral result (solid dots) over the entire range 
of dissipation. The zeroth-order result for (a^) in the the polaron frame (crosses) tends 
to overestimate the correction, whereas the variational frame result (diamonds) provides at 
least a qualitatively correct description. 

B. Slow Bath, ojc < A 

Fig. [2] displays the opposite case, when the bath is slow as compared to the tunneling rate, 
ujc < A. At small and intermediate 7, it can be seen that the polaron method with 2nd-PT 
fails to predict the correct behavior, while the usual 2nd-PT result in the original frame 
agrees well with the exact result. As with the case of the fast bath, at large 7, the polaron 
method provides an accurate description, while the original frame 2nd-PT breaks down. 
The variational polaron method (with 2nd-PT) interpolates between these two methods, 
providing accurate results over a large range of the dissipation strength. The failure of the 
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FIG. 1: (Color online) Fast bath, ujc > A. Comparison of the approximate methods with the 
exact path integral results as a function of the dissipation, 7, for Uc = 5e, and A = 3e. Plotted are 
the values of (a^) from the zeroth order density matrix in the polaron frame (light blue crosses) 
and the variational polaron frame (dark blue diamonds), as well as from the density matrix in the 
original frame with the second order correction (red dashed line), the polaron frame (purple empty 
circles) and the variational polaron frame (solid green line). The exact path integral results are 
shown as filled (black) dots. 

full polaron method is due to the fact that the bath oscillators are sluggish and are not able 
to fully dress the system. Therefore the full polaron displacement is no longer appropriate. 
It can also be seen that the second order correction in the full polaron frame is huge at 
small 7 (the difference between the crosses and open circles) . This should cast doubt on the 
validity of 2nd-PT in this case since the perturbative correction should be small. 

It can also be observed that there is a discontinuity in the variational result (both with 
and without 2nd-PT) at the critical point of 7 10.6. The variational approach exhibits 
a rather abrupt transition from a small transformation [F{(jj) <C 1] to the full polaron 
transformation [F{(jj) ^ 1]. This discontinuity, which is an artifact of the transformation 
rather than a physical phase transition, has been predicted by Silbey and Harris^^'^^ for a 
slow bath. The discontinuity comes from solving the self-consistent equation in Eq. ( IT4I ). 
Over a certain range of dissipation strengths, there exist multiple solutions to the self- 
consistent equation. According to the variational prescription, the solution with the lowest 
free energy is selected. This causes a "jump" in the solution, as depicted in Fig. [31 It is also 
observed in Fig. [2] that the variational polaron result is least accurate around the transition 
point. Therefore, it will be worthwhile to look for a better variational criterion that removes 
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this discontinuity, which can hopefuUy provide a uniformly accurate solution. 
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FIG. 2: (Color online) Slow bath, ujc < A. Comparison of the approximate methods with the exact 
path integral results as a function of the dissipation, 7, for Uc = 1.5e, and A = 5e. Plotted are 
the values of (a^) from the zeroth order density matrix in the polaron frame (light blue crosses) 
and the variational polaron frame (dark blue diamonds), as well as from the density matrix in the 
original frame with the second order correction (red dashed line), the polaron frame (purple empty 
circles) and the variational polaron frame (solid green line). The exact path integral results are 
shown as filled (black) dots. 




B 

FIG. 3: (Color online) Plot of the expression = B — e~^'^ ~ ^^^^ coth(/3a;/2)^ rpj^^ solutions 
to the self consistent equation Eq. (fT4l) are the points when = 0. At 7 = 9.5 there is only one 
solution to the self-consistent equation. At 7 = 10, multiple solutions start to develop, but the 
solution with the lowest free energy (denoted by the empty circle) is chosen. At the critical point 
of 7 = 10.6, there is a "jump" in the lowest free energy solution which causes a discontinuity in 
the transformation. 
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C. Adiabatic Bath, cj^ < /S, A 




FIG. 4: (Color online) ujc <C A. Comparison of the approximate methods with the exact result 
[from Eq. (|26]) ] as a function of the dissipation, 7, for ujc — 0.1, and A = 1. Plotted are the values 
of {dz) from the density matrix with the second order correction in the original frame (dashed red 
line), the polaron frame (empty blue circles) and the variational polaron frame (solid green line). 
The exact results are shown as filled (black) dots. 

In the adiabatic limit {uJc <C /3) the exact solution to the equilibrium state of the system 
can be obtained analytically. In this regime, the partition function is given by — 

Zs = ^^^P{ - ^ + [2cosh[/3V(W + (^ + e/2)2]] }, (26) 

where x is the bath correlation function [in the original frame with F(cl;) = 0] in the adiabatic 
limit, X = Cfz = f^. The expectation value, (a^), can be obtained from the partition 
function via the following relation 

= -||ln2.. (27) 

In the regime where <C A, the transition from F{(jj) = to F{(jj) = 1 in the variational 
method is sharp, as seen in Fig. [H Before the transition, the variational polaron result 
coincides with the exact result and that of perturbation theory in the original frame. The 
full polaron result fails to give the correct results, and even predicts the wrong limiting 
behavior as 7 ^ 0. After the transition, the variational result deviates from the exact result 
and becomes essentially the same as the full polaron result. As 7 increases, results from 
both methods approach the exact result while the untransformed 2nd-PT breaks down as 
seen before. 
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D. Relative Errors 



To get a better perspective of how the accuracy of 2nd-PT in different frames depends 
on the properties of the bath, we calculate the relative errors over the entire range of the 
bath parameters. The relatives errors are defined as 



where the subscripts "Pert" and "PI" denote the perturbative calculation and path integral 
calculation respectively. Fig. [5] displays the respective errors for the three methods as a 
function of the cut-ofl[ frequency and the coupling strength. As seen in Fig. [5](a), the usual 
2nd-PT without transformation breaks down at large 7. It is also less accurate when the 
cut-off frequency is small, which corresponds to a highly non-Markovian bath. On the other 
hand, the 2nd-PT in the full polaron frame fails at small 7 and ujc [see Fig. [5](b)]. These 
two approaches provide complementary behavior as a function of the coupling strength; the 
polaron method is essentially exact for large 7, while the usual 2nd-PT is exact for small 7. 
The variational calculation is valid over a much broader range of parameters [see Fig. [5](c)], 
and essentially combines the regimes of validity of the full polaron result and 2nd-PT in the 
original frame. It is only slightly less accurate in the slow bath regime around the region 
where the discontinuity appears that was discussed above. 

IV. CONCLUSIONS 

In conclusion, we have provided a thorough assessment of the accuracy of the polaron 
and variational polaron methods. We compared the second order perturbation results in 
the polaron and variational polaron transformed frames with numerically exact path inte- 
gral calculations of the equilibrium reduced density matrix. Focusing on the equilibrium 
properties allowed us to systematically explore the whole range of bath parameters without 
making any additional approximations as is generally required to simulate the dynamics. 
As a function of the system-bath coupling, it is found that the standard perturbation result 
without the polaron transformation is accurate for small coupling, while the polaron result 
is accurate in the opposite regime of strong coupling. The variational method is capable of 
interpolating between these two limits. It is valid over a much broader range of parameters 
and is only slightly less accurate around the region where the discontinuity appears. As the 




(28) 
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FIG. 5: (Color online) At A = 3, the relative errors of the second order perturbation theory as 
defined in Eq. (128]) in (a) the original frame, (b) the full polaron frame, and (c) the variational 
polar on frame. 

relaxation time of the bath becomes longer leading to more non-Markovian character, all 
three of the perturbation methods are seen to be less accurate. 
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Appendix A: Imaginary Time Path Integral 

For Hamiltonians such as the spin-boson model where the bath is harmonic, the trace 
over the bath degrees of freedom in Eq. (fT5l) may be performed analytically. In the 
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path integral formulation this procedure leads to the well-known Feynman- Vernon influ- 
ence functional.— Using a Hubbard-Stratonovich transformation,— i^^i^iSl it was shown 
that the influence functional may be unraveled by an auxiliary stochastic fleld. The ensuing 
imaginary time evolution may then be interpreted as one governed by a time-dependent 
Hamiltonian. Explicitly, the spin-boson model may be equivalently expressed as 

HiT) = ^a, + ^a, + a,aT). (Al) 

All of the eflFects of the bath are accounted for by the colored noise term, ^(r), which obeys 
the autocorrelation relation, 

mCir')) = Ci'Xr - r') , (A2) 

where Ci^^(r) is the correlation function given in Eq. I HM with F{(jj) = 0. The trace over the 
bath that was present in the original path integral formulation now corresponds to averaging 
the imaginary time dynamics over realizations of the noise. The auxiliary fleld is simply an 
efl5cient method of sampling the influence functional. 

In practice, a sample of the reduced density matrix is propagated to the imaginary time 
/3, where the time steps, 5r, are determined by 

ps{r + 5T)=exp -5tH{t) ps{r) , (A3) 

with the initial condition, p(0) = /. The primary beneflt of this approach is that it generates 
the entire reduced density matrix from a single Monte Carlo calculation. Additionally, any 
form for the spectral density of the bath, J{uj)^ may be used. In our calculations, 10^ (at 
small 7) to 10^^ (at large 7) Monte Carlo samples are needed to achieve convergence. 
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